function [U_hat,m_hat,f_hat,s_hat,J_hat,W_hat,Gdest_hat,Gunc_hat,meanprices_hat,Jf_hat ,ccp_hat ,lambda_hat, lambda_f_hat] = efficient_equil(s_ineff,f_ineff,U_ineff,J_ineff,Jf_ineff,Gdest_ineff,meanprices_ineff,cobb_par,gamma_new,cs,cu,c_inv,sigma,beta,ksi,g_euler,r,delta,K_obs,E,S_obs,ineff,sigma_exp)
% This function computes the constrained efficient allocation. This is
% implemented by imposing the efficient prices (see eq 31 in the paper)

g0 = 0.5;

g1 = 0.8;

G = 8;

M = size(f_ineff,1);

%
% Initialize starting point
%
f_init = f_ineff;
s_init = s_ineff;
U_init = U_ineff;
J_init = J_ineff;
Gdest_init = Gdest_ineff;
meanprices_init = meanprices_ineff;
Jf_init = Jf_ineff;
share = 1;

%
% Initialize convergence parameters
%

diff_ships = 100;

iter_ships = 1;

eps = exp(-5);

while abs(diff_ships)>0.001 && iter_ships<1000

    
    s_0 = share.*s_init;

    f_0 = f_init;

    U_0 = U_init;

    J_0 = J_init;

    Gdest_0 = Gdest_init;

    meanprices_0 = meanprices_init;

    Jf_0 = Jf_init;

    dif = 100;

    iter = 1;
    
    m_0 = zeros(M,1);
    
    while dif>eps && iter<100000

        %
        % Compute the number of matches, given guesses for the number of
        % searching ships s_0 and the number of searching exporters f_0 
        %
        for i = 1:M

            gamma_i = cobb_par(i,2);

            m_0(i) = cobb_par(i,1);
            m_0(i) = m_0(i).*(s_0(i).^(1-gamma_i));
            m_0(i) = m_0(i).*(f_0(i).^(gamma_i));

            if m_0(i) >min(s_0(i),f_0(i))

                m_0(i) = min(s_0(i),f_0(i)); 
            end

        end

        lambda_f_0 = (m_0./f_0);

        lambda_0 = (m_0./s_0);

        %
        % Compute the ships' value function, given guesses for matching
        % rates and prices
        %

        max_iter_inner = 10000; 

        [ U_1,W,J_1 ] = Vfunc_ship( Gdest_0,cu,cs,ksi,lambda_0,g_euler,meanprices_0,M,max_iter_inner,beta,U_0,sigma );
        ccp_sim = ccp_equil(W,M);

        %
        % Update the guess for the efficient prices
        %
        omega = r - beta.*delta.*Jf_0 + (W - J_1)./sigma;

        meanprices_1 = (J_1 - W)./sigma + (1-gamma_new).*nansum(Gdest_0.*omega,2);

        meanprices_1(eye(size(meanprices_1))==1) = nan;
        
        %
        % Compute exporters' value function and choice probabilities 
        %
        
        Jf_1 = zeros(M,M);

        for i=1:M

            alpha = (1 - delta.*beta.*(1-lambda_f_0(i)));

            Jf_1(i,:) = lambda_f_0(i)*(r(i,:)-meanprices_0(i,:))-c_inv(i);

            Jf_1(i,:) = Jf_1(i,:)./alpha;

        end

        Jf_1(eye(size(Jf_1))==1)=nan;

        vvs = cat(2,zeros(15,1),Jf_1-K_obs)./sigma_exp;
        c = max(vvs,[],2);
        mm = min(vvs,[],2);
        Ma = abs(mm);
        c(Ma>c) = mm(Ma>c);
        g = exp(vvs -c);
        Gunc = g./nansum(g,2);
        if (sum(Gunc(:,1)>0.999)>1)
            vvs = (Jf_1-K_obs)./sigma_exp;
            c = max(vvs,[],2);
            mm = min(vvs,[],2);
            Ma = abs(mm);
            c(Ma>c) = mm(Ma>c);
            g = exp(vvs -c);
            g = g./nansum(g,2);
            Gunc(Gunc(:,1)>0.999,2:16) = 0.001*g(Gunc(:,1)>0.999,:);
            Gunc(Gunc(:,1)>0.999,1) = 0.999;
        end

        Gdest_1 = Gunc(:,2:end)./(1-Gunc(:,1));

        %
        % Update the number of ships and exporters searching
        %
        
        
        s_1 = zeros(M,1);

        for i = 1:M

            s_1(i) = nansum(ccp_sim(:,i).*(s_0-m_0)) + nansum(Gdest_1(:,i).*m_0);

        end

       
        f_1 = delta.*(f_0-m_0) + E.*(1-Gunc(:,1));

        f_1(f_1<0) = 0;
        %
        % Check convergence
        %
        dif_s = (s_1-s_0);

        dif_f = abs((f_0-f_1)./(1+f_0));

        dif_U = abs(U_0-U_1);

        dif_p = abs(meanprices_1(:)-meanprices_0(:))./(1+abs(meanprices_0(:)));

        dif_G = abs(Jf_0(:) - Jf_1(:));

        dif = [dif_s;dif_f;dif_U;dif_p(~isnan(dif_p));dif_G(~isnan(dif_G))];

        dif = max(abs(dif));

   
        coeff = 1-(1./(G+iter^g0));

        f_0 = coeff.*f_0 + (1-coeff).*f_1;

        s_0 = coeff.*s_0 + (1-coeff).*s_1;

        U_0 = coeff.*U_0 + (1-coeff).*U_1;
        
        m_0 = zeros(M,1);
        
        for i = 1:M

            gamma_i = cobb_par(i,2);

            m_0(i) = cobb_par(i,1);
            m_0(i) = m_0(i).*(s_0(i).^(1-gamma_i));
            m_0(i) = m_0(i).*(f_0(i).^(gamma_i));

            if m_0(i) >min(s_0(i),f_0(i))
                m_0(i) = min(s_0(i),f_0(i)); 
            end

        end

        coeff = 1-(1./(G+iter^g1));

        Jf_0 = coeff.*Jf_0 + (1-coeff).*Jf_1;

        Gdest_0 = coeff.*Gdest_0 + (1-coeff).*Gdest_1;
        meanprices_0 = meanprices_1;
        iter = iter+1;

    end
    %
    % Check convergence of the total number of ships in the market
    %
    s_tot = tot_fleet(s_0,m_0,ccp_sim,Gdest_1,ksi);

    diff_ships = (s_tot-S_obs)./S_obs;
    
    if diff_ships>0

        share = share - abs(diff_ships)/(iter_ships^0.5); 

    elseif diff_ships<0

         share = share + abs(diff_ships)/(iter_ships^0.5);
    end

    

    iter_ships = iter_ships+1;
end

%
% Save the equilibrium allocation
%

U_hat = U_0;

m_hat = m_0;

f_hat = f_0;

s_hat = s_0;

J_hat = J_1;

W_hat = W;

Gdest_hat = Gdest_1;

Gunc_hat = Gunc;

meanprices_hat = meanprices_1;

Jf_hat = Jf_0;

ccp_hat = ccp_sim;

lambda_hat = lambda_0;

lambda_f_hat = lambda_f_0;


end

